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1.  Introduction 

Modelling  of  the  inverse  background-error  correlations 
(iBECs)  of  random  fields  by  differential  operators  has  gained 
considerable  attention  in  recent  years,  primarily  due  to 
computational  efficiency  of  their  implementation  in  the 
iterative  minimization  algorithms  used  in  variational  data 
assimilation  (e.g.  Xu,  2005;  Pannekoucke  and  Massart,  2008; 
Mirouze  and  Weaver,  2010).  Of  particular  interest  are  the 
iBEC  models  described  by  positive-definite  polynomials  of 
the  diffusion  operator 

D  =  Vi’V ,  (1) 

where  v  is  the  symmetric  positive-definite  spatially  varying 
diffusion  tensor.  This  type  of  iBEC  model  is  attractive  for 
several  reasons:  (a)  it  allows  straightforward  control  of 
inhomogeneity  and  anisotropy  via  the  diffusion  tensor;  (b) 
it  is  computationally  competitive  in  many  applications; 
and  (c)  it  is  easier  to  develop  with  regard  to  keeping 
the  positive-definiteness  property  of  the  BEC  operator. 
In  contrast,  in  the  traditional  approach  of  the  ‘direct* 


correlation  modelling  where  spatial  correlations  are  specified 
by  prescribed  analytical  functions,  care  should  be  taken  to 
maintain  positive  definiteness  of  the  respective  correlation 
operator,  especially  in  anisotropic  and/or  inhomogeneous 
cases  (e.g.  Gaspari  et  al .,  2006;  Gregorl  et  al. ,  2008). 

Because  of  the  above-mentioned  properties,  polynomials 
of  D  have  been  extensively  used  for  approximating  Gaussian 
shaped  BECs  by  either  explicit  (e.g.  Derberand  Rosati,  1989; 
Egbert  et  al, ,  1994;  Weaver  and  Courtier,  2001;  Weaver 
et  al. ,  2003)  or  implicit  (Ngodock  et  al ,  2000;  DiLorenzo 
et  aly  2007)  integration  schemes.  In  the  latter  case,  the  BEC 
operator  is  obtained  via  iBEC  representation  by  a  binomial 
of  D,  Second-order  polynomials  of  D  were  considered 
recently  by  Hristopulos  (2003)  and  Hristopulos  and  Elogne 
(2007,  2009)  for  construction  of  the  correlation  models 
for  geostatistical  and  other  applications.  Our  research 
(Yaremchuk  et  al .,  2011)  also  indicates  that  low-order 
iBEC  models  can  provide  extra  computational  savings 
in  three-dimensional  variational  (3D-Var)  analysis  while 
keeping  the  predictive  skill  of  oceanographic  assimilation 
systems.  A  comprehensive  treatment  of  representing  the 
iBECs  by  the  polynomials  of  D  was  given  by  Xu 
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(2005),  who  obtained  Taylor  expansions  approximating 
the  propagator  of  the  diffusion  equation  in  one,  two 
and  three  dimensions  and  obtained  recursive  relations 
for  the  polynomial  coefficients  associated  with  spatially 
homogeneous  correlation  functions. 

It  should  be  noted  that,  among  the  numerous  classes 
of  analytic  correlation  functions  (CFs),  only  a  few  can 
be  effectively  implemented  within  the  operator-based 
approach,  because  the  corresponding  iBEC  is  represented 
by  a  polynomial  of  D.  In  many  geophysical  applications, 
however,  details  of  the  shape  of  a  CF  are  poorly  known 
because  of  insufficient  statistics.  As  a  consequence,  in  most 
cases,  heuristic  CFs  can  be  adequately  approximated  using 
only  one  or  two  scalar  parameters  that  define  a  family  of 
analytic  correlation  functions.  Therefore,  analytic  CFs,  with 
inverses  that  can  be  described  by  low-order  polynomials  of 
D,  are  of  significant  practical  interest 

Based  on  the  results  of  the  recent  studies,  this  note  presents 
analytical  expressions  for  the  CFs  corresponding  to  two  types 
of  the  iBEC  models:  the  first  type  is  the  mth-order  binomial 
of  D,  which  approximates  the  Gaussian-shaped  CF;  and  the 
second  type  is  a  quadratic  function  of  D  that  is  capable 
of  reproducing  negative  lobes  in  the  CFs.  The  obtained 
CFs  generalize  earlier  results  of  Hristopulos  and  Elogne 
(2007;  hereinafter  HE07)  and  Mirouze  and  Weaver  (2010), 
and  may  facilitate  practical  design  of  the  cost  functions 
in  variational  data  assimilation  problems,  as  they  give 
explicit  relationships  between  the  shape  of  the  CFs  and  the 
structure  of  the  corresponding  iBEC  operators  in  the  analytic 
form. 

2.  CFs  generated  by  the  polynomials  of  the  homogeneous 
diffusion  operator 

Consider  an  anisotropic  homogeneous  diffusion  operator 
(1)  in  R n,n=  1,...,3  with  x  €  R"  representing  points 
in  the  physical  space  and  k  representing  points  in 
the  wavenumber  space.  By  the  appropriate  coordinate 
transformation  (e.g.  Xu,  2005;  HE07)  the  problem  can 
be  reduced  to  considering  isotropic  BEC  operators  of  the 
form  F(  —  A),  where  F  is  a  positive  function  and  A  is  the 
Laplacian  operator.  In  the  general  case  of  an  inhomogeneous 
operator,  such  transformation  cannot  be  found,  but  local 
transformations  of  this  type  can  be  useful  in  constructing 
the  BEC  operator. 

In  this  note,  the  following  two  classes  of  the  iBEC 
operators  are  considered:  the  first  is  represented  by  the 
binomial 

B-'  =  (l  — aoA)",  (2) 

and  the  other  by  the  second-order  polynomial  in  A: 

B-1  =  I  —  a,A  +a2A2.  (3) 

Here  I  is  the  identity  operator  and  or,  are  the  real  numbers, 
constrained  by  the  positive  definiteness  requirement  of 
B_1.  In  the  binomial  case  (2),  this  constraint  is  cto  >  0. 
For  the  quadratic  polynomial  (3),  the  positive-definiteness 
requirement  can  be  taken  into  account  explicitly  by 
diagonalizing  B~*  via  the  Fourier  transform.  In  the  Fourier 
representation,  B“ 1  acts  as  multiplication  by  the  polynomial 
inJt2  =  |  k|2,  and  the  positive-definiteness  property  translates 
into  the  requirement  for  the  spectral  polynomial 

B-,(Jt)=  1  +a,*2+a2Jt<  (4) 


to  be  positive  for  all  k2  =  |fc|2  (e.g.  Reed  and  Simon,  1975). 
This  constraint  is  equivalent  to  the  statement  that  the 
polynomial  in  the  right-hand  side  of  (4)  must  not  have 
real  positive  roots.  Since  we  are  considering  biquadratic 
polynomials  with  real  coefficients,  these  roots  are  symmetric 
with  respect  to  both  real  and  imaginary  axes.  Thus,  without 
loss  of  generality  (except  for  the  special  case  of  imaginary 
roots,  which  is  treated  later),  B~l  (k)  can  also  be  represented 
in  the  form 

B~'(k)  =  y{a2  +  (Jt  -  fo)2}{«2  +  (Jfc  +  b )2} ,  (5) 

where  a  and  b  are  real  numbers  defining  the  inverse 
decorrelation  scales  of  the  covariance  operator,  and  y  = 
(a2  +  fr2)-2.  The  correspondence  between  a\ya2  and  ayb 
can  easily  be  established: 

cr  j  =  2(a2  —  b2)(a2  +  fr2)~2;  otj  =  (u2  +  b2)~2  (6) 

Compared  to  the  spectral  representation  (4)  considered  in 
HE07,  the  representation  (5)  has  the  advantage  that  its  free 
parameters  are  not  constrained  by  the  positive-definiteness 
requirement.  The  reciprocal  of  B~*(k)  provides  the  spectral 
representation  of  the  BEC  operator: 

B(k )  =  [y (a2  +  {k-  />)2}{«2  +  (*  +  fc)2}]"1 .  (7) 

In  a  special  case  when  both  roots  are  on  the  imaginary  axis, 
the  diagonal  of  B  1  can  be  represented  by 

B"l(k)«y(<i2  +  it2)(fc2+  it2).  (8) 

where  y  =  (ab)~2  and  the  weighting  factors  before  the 
Laplacians  are  given  by 

c?i  =(a2  +  b2)(ab)~2y  ai  *=(ab)~2.  (9) 

The  corresponding  spectrum  can  be  reduced  to  the 
difference  of  the  respective  first-order  (one-parameter) 
spectra 


y(a2  +  k2)^  +  k2) 


b2  -a2  \a2  +  k2  P  +  k2/ 

considered  in  the  next  section. 

Because  of  homogeneity,  the  matrix  elements  of  B  depend 
only  on  the  distance  r  =  |x|  from  the  diagonal.  They  can 
be  computed  by  applying  the  inverse  Fourier  transform  to 
B  (k): 

B"(r)  =  J  B(k)  exp(— lJfcx)  dJt .  (10) 

R" 

By  integrating  over  the  directions  in  R”  (see  the  Appendix), 
(10)  can  be  reduced  to 

oo 

Bn(r)  =  (27T)-n/1  J  B(Jt))t"- 1  (Jtr)-|/,( Jtr)  dJt ,  (11) 

0 

where  J  denotes  the  Bessel  function  of  the  first  kind 
and  s  =  n/2  —  1.  The  respective  matrix  elements  of  the 
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correlation  operator  (correlation  functions)  are  obtained  by  In  the  limiting  case  m  — ►  oo,  the  correlation  functions  (15) 
normalization:  take  the  Gaussian  form: 


C"(r)  =  B"(r)/BM(0).  (12) 

In  practical  applications,  the  diffusion  operator  is  not 
homogeneous,  and  the  analytic  representations  (4)-(ll) 
cannot  be  obtained.  However,  the  action  of  the  BEC  operator 
on  a  state  vector  can  be  computed  numerically  at  a  relatively 
low  cost.  The  major  problem  with  such  modelling  is  the 
efficient  estimation  of  the  diagonal  elements 

B"(x,x)  =  j  Bn(x,y)S(x-y)dy,  (13) 

R" 

which  are  necessary  to  rescale  B  to  have  its  diagonal  elements 
equal  to  unity.  In  practice,  the  rescaling  factors  Nn(x)  are 
defined  as  reciprocals  of  Btt(x,  x). 

Taking  the  integral  in  (13)  numerically  is  expensive, 
because  the  convolutions  with  the  8 -functions  have  to 
be  performed  at  all  numerical  grid  points  x.  However, 
reasonable  approximations  for  Nn(x)  can  be  obtained  by 
using  the  homogeneous  analytical  versions  of  (13)  (e.g. 
Purser  et  «/.,  2003;  Mirouze  and  Weaver,  2010).  Therefore, 
analytical  formulae  describing  homogeneous  BEC  operators 
are  of  significant  practical  interest.  Another  benefit  of  the 
analytical  models,  is  their  ability  to  provide  guidance  in  the 
design  of  the  correlation  functions.  In  the  case  considered, 
the  type  of  spectral  polynomial  defines  the  CF’s  shape 
as  a  function  of  a,.  Conversely,  it  provides  the  values  of 
a,  after  the  CF  parameters  are  (optimally)  fitted  to  the 
available  data.  In  this  note,  two  types  of  such  polynomials  are 
considered:  the  first  type  describes  power  approximations 
of  the  Gaussian-shaped  CF,  and  the  second  is  a  general 
second-order  polynomial  given  by  (5)  and  (8). 

3.  Power  approximations  of  the  Gaussian-shaped  CF 

An  important  family  of  one-parameter  correlation  spectra 
provides  approximations  to  the  Gaussian-shaped  correlation 
function: 

<"> 


C”oc  =  exp(-r2/2a2);  n  =  1,...  (17) 

Consecutive  approximations  of  the  Gaussian  CF  by  (15)  are 
shown  in  Figure  1.  It  is  remarkable  that,  when  m  =  1,  the 
correlation  functions  (15)  have  singularities  at  p  =  0  in  both 
two  and  three  dimensions  (also  Table  I).  This  means  that  in 
the  continuous  case  the  first-order  approximations  become 
invalid  when  n  >  1.  Numerically,  however,  the  correlation 
functions  do  exist  for  n  >  1,  but  their  decorrelation  scale  is 
limited  by  the  grid  size  <5  (the  corresponding  CF  is  shown 
by  the  dotted  line  in  Figure  1(a)).  This  occurs  because  the 
numerical  analogue  of  the  <5-function  is  never  singular,  but 
has  a  finite  amplitude  inversely  proportional  to  the  volume 
of  a  grid  cell,  therefore  resulting  in  a  finite  value  of  the 
convolution  (13)  even  if  it  is  infinite  in  the  continuous  case. 
After  normalization  by  that  finite  value,  the  CF  is  1  at  r  =  0, 
but  its  effective  decorrelation  scale  remains  proportional  to 
the  local  grid  size  8  if  a  >$>  8. 

It  is  also  noteworthy  that  the  rrtth-order  correlation 
functions  in  3D  coincide  with  the  l)th-order  CFs 
in  ID.  In  particular,  the  ID  second-order  autoregression, 
or  SOAR  function,  widely  used  in  operational  analyses, 
corresponds  to  the  third-order  approximation  of  the 
Gaussian  function  in  3D. 

Figure  1(a)  shows  that  low-order  power  approximations 
(14)  underestimate  the  decorrelation  scale  a  of  the  target 
Gaussian  function.  This  unpleasant  property  can  be 
corrected  by  optimizing  the  value  of  a  in  (14)  to  obtain 
the  best  fit  with  the  Gaussian  CF.  Because  the  Gaussian 
and  its  approximating  functions  are  both  positive  and  have 
similar  shapes,  a  reasonable  optimization  criterion  is  to  set 
their  integral  decorrelation  scales  equal  to  each  other: 

00  oo 

/  C”Jp)dr=-^LfcXy)dy 

0  °  (18) 

oo  _  v  ' 

=/”,,(“^)  t‘Vf- 

0 


Although  the  binomial  approximation  (14)  converges  fast 
enough  (e.g.  Abramowitz  and  Stegun,  1972),  only  small 
values  of  m  are  of  practical  interest.  In  this  section  we  derive 
the  binomial-generated  CFs  and  the  correction  coefficient 
needed  for  efficient  approximation  of  the  Gaussian  CF  when 
m  is  small. 

Substituting  (14)  into  (11),  integrating  over  k  and 
normalizing  the  result  by  BM(0)  yields  the  correlation 
functions  of  the  Matern  family  (Stein,  1999)  enumerated 
by  5  =  m  —  n/2  and  scaled  by  a*  =  a/y/2m: 


CnJp)  = 


p%(p) 

2!-T(s)’ 


(15) 


where  p  =  r/a*,  T  is  the  gamma  function  and  K  stands 
for  the  modified  Bessel  function  of  the  second  kind 
(e.g.  Abramowitz  and  Stegun,  1972).  The  respective 
normalization  factors  are 

K  =  ^(2^.)".  (16) 


Expression  (18)  shows  that  aopt  =  f”a,  with  the  rescaling 
coefficient 


r(s) 

r(s+ 1/2) 


yfm. 


(19) 


The  values  of  ££  for  me  4  and  their  respective 
approximation  errors  are  assembled  in  Table  1. 

The  coefficients  £”  along  with  relationship  (14)  provide 
an  expression  for  estimating  a0  in  the  binomial  iBEC  model 
(2)  which  approximates  the  Gaussian-shaped  correlation 
function  with  a  given  radius  a: 

<*o  =  (t>)2/2m  (20) 


4.  Two-parameter  correlation  functions 


In  the  general  case  of  the  two-parameter  approximation 
(5)  there  are  two  complex  roots  located  symmetrically  with 
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fa)  (b) 


(c) 


Figure  1.  (a)  Power  approximations  (14)  of  the  Gaussian  CF  in  two  dimensions  (r*  =  2).  The  CF  for  m  =  1  is  shown  by  the  dott«i  line  for  tht  numerical 
realization  with  the  grid  step  &  =  a/4.  (b)  Same  approximations,  but  with  optimally  adjusted  correlation  radii  for  various  combinations  of  m  and  n.  (c) 
Differences  between  the  Gaussian  CF  and  its  approximation  shown  in  (b).  The  horizontal  axes  are  scaled  by  a. 


Table  I.  Correlation  functions  associated  with  the  power  approximations  (14)  of  the  Gaussian  CF  in  n  dimensions. 
The  CFs  for  n  =  1,3  are  rewritten  in  terms  of  elementary  functions  for  convenience.  The  correlation  radius  adjustment 
coefficients  £  ”  are  shown  below  the  formulae  together  with  (in  bold)  the  corresponding  relative  errors  in  approximation 

of  the  Gaussian  CF. 


n  =  1 

n  =  2 

n  =  3 

m  =  1 

exp(-p) 

Ko(p) 

exp  (— p)/p 

yfn  0.33 

— 

— 

m  =  2 

(1  +  p)exp(-p) 

pKi(p) 

exp(-p) 

\fnji  0.13 

y/Sjn  0.19 

v/2tt  0.33 

m  =  3 

(1  +  P  +  P2/3)  exp(— p) 

P2K2(p)/ 2 

(1  +p)exp(-p) 

Vrin/S  0.08 

7I673F  0.10 

v/37t74  0.13 

m  =  oo 

exp(— r2/2«2) 

exp(— r2/2o2) 

exp(— i^/la1) 

which  isn/2b  away  from  the  origin  for  n  =  3,  and  depends 
on  both  a  and  b  for  n  =  1 .  Although  analytical  expressions 
are  quite  different  for  n  =  2  and  n  =  3,  the  behaviour  of  the 
CFs  is  rather  similar.  In  the  1 D  case,  the  first  zero  is  somewhat 
farther  away  from  the  origin  and  the  CF  is  less  damped. 

In  the  special  case  (8),  the  CFs  can  be  expressed  via  the 
differences  of  the  first-order  CFs  (15)  discussed  in  section  3: 


respect  to  imaginary  axis.  Plugging  (7)  into  (11),  integrating 
over  ky  and  renormalizing  yields  the  following  CFs: 

C"(o,  b,  r)  =  2-  {z%(z)-z%(z)} ,  (21) 

MM 

where  z  =  (a- b  ib)ry  the  overline  denotes  the  complex 
conjugation,  5=1—  n/ 2,  and  the  coefficients  fin  are 


rjfp  -  a  =  2arctan(:)  •  (22) 

Note  that,  despite  a  seemingly  complex-valued  expression  in 
the  right-hand  side  of  (21),  its  imaginary  part  is  identically 
zero.  Similar  to  the  binomial  case,  1 D  and  3D  two-parameter 
CFs  can  also  be  expressed  in  terms  of  elementary  functions: 


C\aybyr)  = 


C3(a,b,  r)  =  exp(— ar) 


\/a2  T  b2 

b 

x  exp(—  ar)  cosjfcr— arctan  , 

sin  (br) 


hr 


(23) 

(24) 


Figure  2  shows  the  dependence  of  the  correlation 
functions  (21)  on  the  magnitude  of  b  for  various  n.  As 
can  be  seen  from  (23)  and  (24),  CFs  in  ID  and  3D  have 
equidistant  zeros  separated  by  nib  except  for  the  first  zero 


Cl(iiyby  r)  = 
C2(Uyby  r)  - 

C3(«,  byr)  = 


aexp(-br)  —bexp(-ar) 
a  —  b 

K0(ar)  -  K0(br) 
log (b/a) 

exp(— ar)  —  exp(— br) 
(b-a)r 


(25) 

(26) 
(27) 


In  the  expressions  (25),  (27)  the  representations  of  the 
Bessel  functioasin  terms  of  elementary  functions  were  used. 
Also  note  that  C2,3(r)  are  non-singular  at  r  =  0  because  the 
singularities  are  cancelled  out  by  taking  the  difference  in 
the  numerator.  In  this  special  case,  the  second  parameter 
gives  little  extra  freedom  in  adjusting  the  shape  of  the  CFs, 
because  the  resulting  curves  remain  positive  functions  of  r. 
The  extra  degree  of  freedom  can  be  used  to  partly  control, 
for  example,  the  CF  derivative  at  r  =  0. 
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Correlation  Functions  Associated  with  Differential  Operators 


1931 


(a)  (b) 


(c) 


Figure  2.  Two-parameter  CFs  corresponding  to  the  quadratic  inverse  BEC  (5)  with  a  »  I.  The  horizontal  axis  is  scaled  by  a.  Dotted  lines  show  CFs 
corresponding  to  the  case  (8)  with  two  imaginary  roots  of  the  spectral  polynomial. 


The  normalization  constants  for  the  functions  (21)  and 
(25)— (27)  are  respectively 


N]  — 


4  a 

a2  +  b2* 


N2  = 


87 rab 

p2(a2  +  b2)2 


N3  = 


8  Tza 

(a2  +  b2)2 
(28) 


and 

~  ,  2 (a  +  b)  ~2  2n(a2  -  b2)  ~3  4n(a  +  b) 

~ab’  ~  a2b2\og(a/b)'  aW  ' 

(29) 


Equations  (21)  and  (25)-(27)  provide  explicit  expressions 
for  the  CFs  of  the  two-parameter  BEC  model.  Combining 
them  with  the  relationships  (6)-(9)  allows  the  parameters 
of  inverse  BEC  operator  (3)  to  he  computed  after  the  values 
of  a  and  b  are  adjusted  to  experimental  data  using  (21)  or 
(25)— (27). 


5.  Summary  and  discussion 

BEC  modelling  with  diffusion  operators  is  an  efficient 
and  flexible  tool  for  evaluating  matrix-vector  products  of 
large  dimension  which  emerge  in  minimization  algorithms 
of  variational  data  assimilation.  This  note  has  discussed 
analytic  relationships  between  the  parameters  controlling 
the  shape  of  correlation  functions  and  the  polynomial 
coefficients  characterizing  the  structure  of  the  respective 
inverse  BEC  operator.  The  results  may  he  helpful  in 
designing  the  BEC  operators  in  variational  data  assimilation 
algorithms. 

Although  only  homogeneous  operators  in  boundaryless 
domains  were  considered,  the  relationships  ( 15)— ( 16), 
(20)— (29)  may  provide  reasonable  guidance  to  constructing 
more  realistic  BEC  operators,  especially  in  cases  when  the 
typical  scale  of  variability  of  the  diffusion  tensor  is  much 
larger  than  the  local  decorrelation  scale  pc  and/or  most  of  the 
observations  are  separated  from  the  boundaries  by  distances, 
exceeding  pc.  In  a  similar  way,  weak  inhomogeneity  can  be 
introduced  by  variable  scaling  factors  a(x),  b{x ),  and  the 
local  CF  shapes  can  be  assessed  using  (21 )— (27). 

Although  generalizations  of  (21)  for  higher-order  poly¬ 
nomials  are  possible,  this  study  has  been  limited  to  quadratic 
polynomials  for  two  reasons.  First,  the  BEC  operators 
encountered  in  geophysical  fluid  dynamics  applications  are 
rarely  homogeneous  and  observational  statistics  are  usually 
insufficient  to  capture  the  spatial  dependence  of  the  BEC 


structure.  Therefore  experimental  estimates  of  the  BECs  are 
either  limited  to  low-rank  ensemble  estimates  or  have  to  rely 
on  the  very  rough  assumption  of  homogeneity.  Needless  to 
say,  in  the  latter  case  the  structure  of  a  sample  CF  should  be 
elaborated  with  sufficiently  low  detalization  which  can  be 
well  accounted  for  by  a  two- parameter  CF  family.  The  second 
reason  is  that  the  use  of  higher-order  polynomials  consider¬ 
ably  degrades  the  conditioning  of  the  linear  systems  that  are 
being  solved  in  the  assimilation  process  (Yaremchuk  et  ai> 
2011)  and,  therefore,  requires  sophisticated  preconditioners. 

It  should  be  noted  that  similar  problems  have  been 
recently  studied  by  many  authors  (e.g.  Xu,  2005;  Hristopulos 
and  Elogne,  2007,  2009;  Mirouze  and  Weaver,  2010).  In 
particular,  analytic  formulae  analogous  to  (23),  (24),  (25) 
and  (27)  were  derived  in  a  somewhat  different  setting  by 
HE07  who  considered  iBECs  of  similar  structure.  Xu  (2005) 
analyzed  Taylor  expansions  of  the  Gaussian  BEC  operator 
and  obtained  recursive  relations  for  the  polynomial 
coefficients  associated  with  an  arbitrary  CF.  Mirouze  and 
Weaver  (2010)  also  demonstrated  a  possibility  to  generate 
oscillating  CFs  using  higher-order  polynomials  in  ID. 

The  objective  of  this  note  was  to  present  the  accumulated 
information  in  concise  form  and  to  provide  explicit 
relationships  between  the  polynomial  coefficients  of  the 
iBEC  operators  and  the  corresponding  CF  parameters  that 
can  be  derived  from  experimental  data.  In  addition  to 
this,  coefficients  for  the  power  approximations  of  the 
Gaussian  BEC  operator,  and  the  analytic  expression  (21) 
for  the  two-parameter  model  in  arbitrary  dimension,  have 
been  obtained.  The  latter  includes,  in  particular,  the  2D  case 
formulae  (26),  (28),  (29)  absent  in  HE07,  who  considered 
only  ID  and  3D  cases. 

We  believe  this  note  may  facilitate  further  development  of 
the  BEC  models  in  variational  data  assimilation.  Moreover, 
since  the  described  methodology  can  be  used  for  the 
approximation  of  arbitrary  self-adjoint  operators  with 
positive  spectrum,  results  may  also  find  applications  beyond 
the  BEC  modelling  in  geophysical  inverse  problems. 

Appendix 

I^t  0  be  the  angle  between  x  and  k  in  R"  and  n  >  2.  Then 
the  integral  ( 10)  can  be  rewritten  in  spherical  coordinates  as 

oc 

BT\r)  =  (2n)~nj  B(k)Jexp(  — ikrcosO)^'1  dfcdQ„_i, 
o 

(A.l) 
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where  d£2„_i  is  the  element  of  the  surface  area  of  the 
unit  sphere.  Since  cos#  changes  symmetrically  within  the 
limits  of  integration,  the  imaginary  part  of  the  exponent 
vanishes.  Furthermore,  using  the  identity  d£2„_t  =  dft„_2  • 
sin"-2  6  dfl,  the  integral  (A.l)  can  be  rewritten  as 


oo 

B"{r)  =  (2n)~"jB(k)k"~l  dk  J 

0 

7T 

f° 


2 


x  I  cos(krcosO)  sin”  z9d9. 


(A.2) 


Integration  over  9  (3.715.21  of  Gradshteyn  and  Ryzhik, 
1980)  and  substitution  of  the  formula  for  the  surface  of 
(ti  —  2)-dimensional  unit  sphere  into  (A.2)  yields  (11). 

The  general  relationship  (11)  also  holds  for  n=  1,2 
although  these  cases  require  a  special  (less  complicated) 
treatment. 
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